home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
CD ROM Paradise Collection 4
/
CD ROM Paradise Collection 4 1995 Nov.iso
/
os2
/
adaptor.zip
/
ADAPT.ZIP
/
adaptor
/
examples
/
pde
/
expl21.f
next >
Wrap
Text File
|
1993-06-28
|
24KB
|
714 lines
PROGRAM EXPL2
C***********************************************************************
C *
C AUTHOR: YVONNE LUH, GMD, I1.HR *
C *
C POSTFACH 1319 *
C W-5205 ST. AUGUSTIN *
C *
C***********************************************************************
C***********************************************************************
C fortran 90 - Version
C 28/4/92
C.......................................................................
C EXPLICIT ALGORITHM: STANDARD 2ND ORDER DISCRETIZATION
C WAVE EQUATION (2 SPACE DIMENSIONS)
C
C INITIAL BOUNDARY VALUE PROBLEM
C DDU/DDT=DDU/DDX+DDU/DDY
C U(X,Y,0) = F1(X,Y)
C DU/DT(X,Y,O) = F2(X,Y)
C U(X,Y,T) = G(X,Y,T) (X,Y) AT BOUNDARY OF (0,1)X(0,1)
C.......................................................................
C YOU CAN TRY DIFFERENT K,H
C K STEP IN T-DIRECTION
C H STEP IN X- AND Y-DIRECTION
C.......................................................................
C DIFFERENT EXAMPLES CAN BE CALCULATED
C COMMON /NEXAMP /NEXAMP
C 1 DSIN(PI*X)*DSIN(PI*Y)
C * DCOS(DSQRT(2.0)*PI*T)
C * 1000.0D0
C 2 DSIN(20.0*PI*X)*DSIN(20.0*PI*Y)
C * DCOS(20.0D0*DSQRT(2.0D0)*PI*T)
C 3 DEXP(X+Y)*DEXP(DSQRT(2.0D0)*T)
C * 100000.0D0
C 4 DEXP( 3.0D0*(X+Y))
C * DEXP( 3.0D0*DSQRT(2.0D0)*T)
C 5 X**2 + Y**2 + 2.0D0*T**2
C 6 0.0
C DIFFERENT KINDS OF ERRORS CAN BE CALCULATED
C COMMON /NRELER /NRELER
C 1 RELATIVE ERROR =MAX((SOL - U)/SOL)
C 2 ABSOLUTE ERROR =MAX (SOL - U)
C***********************************************************************
C***********************************************************************
c
double precision u0(:,:), u1(:,:), u2(:,:), tx(:,:), ty(:,:)
double precision hu1 (:,:)
double precision hu2 (:,:)
DOUBLE PRECISION TEND,T,K,H,P
DOUBLE PRECISION mxresul, l2resul
INTEGER NP,N
INTEGER NEXAMP,NRELER
INTEGER I, J
c
C.......................................................................
c
C COMMON BLOCKS
C NEXAMP NUMBER OF EXAMPLE
C NRELER RELATIVE OR ABSOLUTE ERROR IS CALCULATED
c
COMMON /NEXAMP /NEXAMP
COMMON /NRELER /NRELER
c
C***********************************************************************
C***********************************************************************
C INPUT OF PROBLEM PARAMETERS
C***********************************************************************
C***********************************************************************
C
READ(5,*) K
READ(5,*) NP
READ(5,*) TEND
READ(5,*) NEXAMP
READ(5,*) NRELER
C
C***********************************************************************
C***********************************************************************
C declaration of allocatable arrays
C***********************************************************************
C***********************************************************************
C
allocate (u0(0:np,0:np), u1(0:np,0:np), u2(0:np,0:np),
+ tx(0:np,0:np), ty(0:np,0:np))
allocate (hu1(0:np,0:np))
allocate (hu2(0:np,0:np))
C
C***********************************************************************
C***********************************************************************
C COMPUTATION OF PARAMETERS
C***********************************************************************
C***********************************************************************
C
N=NP-1
H=1.0D0/DBLE(NP)
P=K/H
C
C***********************************************************************
C***********************************************************************
C calculation of arrays tx and ty
C***********************************************************************
C***********************************************************************
C
!HPF$ INDEPENDENT, LOCAL_ACCESS
do j = 0, np
!HPF$ INDEPENDENT, LOCAL_ACCESS
do i = 0, np
tx(i,j) = dble (i)
ty(i,j) = dble (j)
end do
end do
tx = tx * h
ty = ty * h
c
C***********************************************************************
C***********************************************************************
C OUTPUT OF PARAMETERS
C***********************************************************************
C***********************************************************************
C
WRITE(6,*)'*****************************************************'
WRITE(6,*)'EXPLICIT DISCRETIZATION OF WAVE EQUATION (2D)'
WRITE(6,*)'*****************************************************'
WRITE(6,*) ' COMPUTED EXAMPLE: '
IF (NEXAMP.EQ.1) THEN
WRITE(6,*) ' SIN(PI*X)*SIN(PI*Y)*COS(1.414*PI*T) *1000'
ELSE IF (NEXAMP.EQ.2) THEN
WRITE(6,*) ' SIN(20*PI*X)*SIN(20*PI*Y)*COS(20*1.414*PI*T)'
ELSE IF (NEXAMP.EQ.3) THEN
WRITE(6,*) ' EXP(X+Y)*EXP(1.414*T) * 100000'
ELSE IF (NEXAMP.EQ.4) THEN
WRITE(6,*) ' EXP( 3*(X+Y))*EXP( 3*1.414*T)'
ELSE IF (NEXAMP.EQ.5) THEN
WRITE(6,*) ' X**2 + Y**2 +2*T**2'
ELSE IF (NEXAMP.EQ.6) THEN
WRITE(6,*) ' 0.0 '
ENDIF
IF (Nreler.EQ.1) THEN
WRITE(6,*) ' relative error is calculated'
ELSE IF (Nreler.EQ.2) THEN
WRITE(6,*) ' absolute error is calculated'
ENDIF
WRITE(6,*) '****************************************************'
WRITE(6,9000) ' TIME -STEP K = ', K
WRITE(6,9000) ' SPACE-STEP H = ', H
WRITE(6,9000) ' RATIO P=K/H = ', P
C
IF (P*P.GT.0.5D0) THEN
WRITE(6,*) ' '
WRITE(6,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
WRITE(6,*) ' CFL CONDITION (P*P =< 0.5) NOT FULFILLED'
WRITE(6,*) ' THIS WILL CAUSE LARGE ERRORS IN COMPUTED SOLUTION'
WRITE(6,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
WRITE(6,*) ' '
ENDIF
C
WRITE(6,9001) ' NUMBER OF GRID POINTS PER LINE = ',NP
WRITE(6,9000) ' MAX. T = ', TEND
9000 FORMAT(1X,A17,F20.10)
9001 FORMAT(1X,A35,I5)
WRITE(6,*) '*****************************************************'
C
C***********************************************************************
C***********************************************************************
C COMPUTATION OF INITIAL VALUES IN FIRST TWO TIME STEPS
C***********************************************************************
C***********************************************************************
C
call f1(u0,tx,ty,np)
call sol(u1,tx,ty,np,k)
T=K
C
C***********************************************************************
C***********************************************************************
C TEST OF INITIAL VALUES AT TIME T=K
C***********************************************************************
C***********************************************************************
C
call difmx(u1,tx,ty,np,k,mxresul)
call difl2(u1,tx,ty,np,k,l2resul)
WRITE(6,9020) ' T = ', ' difmx ', ' difl2 '
WRITE(6,9021) t, mxresul, l2resul
9020 FORMAT (1X,A8,2x,a10,2x,a10)
9021 FORMAT (1X,F8.4,2x,D10.4,2x,D10.4)
C
T=K+K
C
C***********************************************************************
C***********************************************************************
C
40 CONTINUE
C
C***********************************************************************
C***********************************************************************
C COMPUTATION OF SOLUTION
C***********************************************************************
C***********************************************************************
C
C IN INTERIOR
C
C overlap U1(:,1:1) with HU1
C
C n = np - 1
hu1(0:np,1:n) = u1(0:np,2:n+1)
hu2(0:np,1:n) = u1(0:np,0:n-1)
!HPF$ INDEPENDENT, LOCAL_ACCESS
do j=1,n
!HPF$ INDEPENDENT, LOCAL_ACCESS
do i=1,n
U2(i,j) = P*P* ( U1(i+1,j) + U1(i-1,j)
+ + Hu1(i,j) + HU2(i,j)
+ - 4.0D0 * U1(i,j) )
+ + 2.0D0 * U1(i,j) - U0(i,j)
end do
end do
C
C.......................................................................
C
C ON BOUNDARY
C
call g(u2,tx,ty,np,t)
C
C***********************************************************************
C***********************************************************************
C COMPUTATION AND OUTPUT OF ERROR OF THE COMPUTED SOLUTION
C***********************************************************************
C***********************************************************************
C
call difmx(u2,tx,ty,np,t,mxresul)
call difl2(u2,tx,ty,np,t,l2resul)
WRITE(6,9021) t, mxresul, l2resul
C
C***********************************************************************
C***********************************************************************
C STORAGE OF NEW SOLUTION
C***********************************************************************
C***********************************************************************
C
U0 = U1
U1 = U2
C
C***********************************************************************
C***********************************************************************
C NEXT TIME STEP
C***********************************************************************
C***********************************************************************
C
T=T+K
IF (T.LT.TEND) GOTO 40
C
C***********************************************************************
C***********************************************************************
C deallocation
C***********************************************************************
C***********************************************************************
c
deallocate (hu2, hu1)
deallocate (ty, tx, u2, u1, u0)
c
C***********************************************************************
C***********************************************************************
C
END
C
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C
C
subroutine DIFMX(U,tx,ty,Np,T,mxresul)
C
C COMPUTES THE MAXIMUM NORM OF THE DIFFERENCE BETWEEN
C SOL (=SOLUTION OF THE BVP) AND THE VALUES IN U
C
C DEPENDING ON NRELER THE RELATIVE ERROR (NRELER=1)
C OR THE ABSOLUTE ERROR (NRELER=2) IS COMPUTED
C
integer np, n
DOUBLE PRECISION U(0:np,0:np), tx(0:np,0:np), ty(0:np,0:np)
DOUBLE PRECISION mxresul
DOUBLE PRECISION uhelp (0:np,0:np)
DOUBLE PRECISION uhelp1 (0:np,0:np)
INTEGER t
INTEGER NRELER
COMMON /NRELER /NRELER
C
N = Np-1
C
call sol(uhelp,tx,ty,np,t)
c
uhelp1 = uhelp
if (nreler.eq.1) then
where ( uhelp(1:n,1:n) .gt. 1.0d-10 )
uhelp1(1:n,1:n) = abs( (uhelp(1:n,1:n)-u(1:n,1:n))
+ / uhelp(1:n,1:n) )
elsewhere
uhelp1(1:n,1:n) = abs( uhelp(1:n,1:n) - u(1:n,1:n) )
endwhere
else if (nreler.eq.2) then
uhelp1(1:n,1:n) = abs( uhelp(1:n,1:n) - u(1:n,1:n) )
endif
c
mxresul = maxval (uhelp1(1:n,1:n))
END
C
C
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C
C
subroutine DIFl2(U,tx,ty,Np,T,l2resul)
C
C COMPUTES THE L2 NORM OF THE DIFFERENCE BETWEEN
C SOL (=SOLUTION OF THE BVP) AND THE VALUES IN U
C
C DEPENDING ON NRELER THE RELATIVE ERROR (NRELER=1)
C OR THE ABSOLUTE ERROR (NRELER=2) IS COMPUTED
C
integer np, n
DOUBLE PRECISION U(0:np,0:np), tx(0:np,0:np), ty(0:np,0:np)
DOUBLE PRECISION l2resul
DOUBLE PRECISION uhelp(0:np,0:np)
DOUBLE PRECISION uhelp1(0:np,0:np)
INTEGER t
INTEGER NRELER
COMMON /NRELER /NRELER
C
N = Np-1
C
call sol(uhelp,tx,ty,np,t)
c
uhelp1 = uhelp
if (nreler.eq.1) then
where ( uhelp(1:n,1:n) .gt. 1.0d-10 )
uhelp1(1:n,1:n) = abs( (uhelp(1:n,1:n)-u(1:n,1:n))
+ / uhelp(1:n,1:n) )
elsewhere
uhelp1(1:n,1:n) = abs( uhelp(1:n,1:n) - u(1:n,1:n) )
endwhere
else if (nreler.eq.2) then
uhelp1(1:n,1:n) = abs( uhelp(1:n,1:n) - u(1:n,1:n) )
endif
c
uhelp1(1:n,1:n) = uhelp1(1:n,1:n) ** 2
l2resul = sum ( uhelp1(1:n,1:n) )
c
l2resul = dsqrt(l2resul) / dble(n)
c
END
C
C
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C
C
subroutine SOL(u,tx,ty,np,T)
C
C**********************************************************************
C
C CONTINUOUS SOLUTION OF THE DIFFERENTIAL PROBLEM
C
C**********************************************************************
C
integer np
DOUBLE PRECISION u(0:np,0:np), tx(0:np,0:np), ty(0:np,0:np)
DOUBLE PRECISION t
DOUBLE PRECISION PI
INTEGER NEXAMP
COMMON /NEXAMP /NEXAMP
PI = DACOS(-1.0D0)
C
C......................................................................
C
IF (NEXAMP .eq. 1) GOTO 1
IF (NEXAMP .eq. 2) GOTO 2
IF (NEXAMP .eq. 3) GOTO 3
IF (NEXAMP .eq. 4) GOTO 4
IF (NEXAMP .eq. 5) GOTO 5
IF (NEXAMP .eq. 6) GOTO 6
C
C......................................................................
C
C NON OSCILLATING SOLUTION
C
1 u = dsin(pi*tx)*dsin(pi*ty)*dcos(dsqrt(2.0d0)*pi*t) * 1000.0d0
GOTO 10
C
C......................................................................
C
C HIGHLY OSCILLATING SOLUTION
C
2 u = dsin(20.0d0*pi*tx)*dsin(20.0d0*pi*ty)
+ * dcos(20.0d0*dsqrt(2.0d0)*pi*t)
GOTO 10
C
C......................................................................
C
C SOLUTION INCREASING SLOWLY IN TIME
C
3 u = dexp(tx+ty) * dexp(dsqrt(2.0d0)*t) * 100000.0d0
GOTO 10
C
C......................................................................
C
C SOLUTION INCREASING RAPIDLY IN TIME
C
4 u = dexp(3.0d0*(tx+ty)) * dexp(3.0d0*dsqrt(2.0d0)*t)
GOTO 10
C
C......................................................................
C
C POLYNOMIAL IN X AND Y AS SOLUTION
C
5 u = tx**2 + ty**2 + 2.0d0*t**2
GOTO 10
C
C......................................................................
C
C ZERO SOLUTION
C
6 u = 0.0D0
GOTO 10
C
C......................................................................
C
10 CONTINUE
END
C
C
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C
C
subroutine F1(u,tx,ty,np)
C
C**********************************************************************
C
C INITIAL CONDITION U(X,Y,0) = F1(X,Y)
C
C**********************************************************************
C
integer np
DOUBLE PRECISION u(0:np,0:np), tx(0:np,0:np), ty(0:np,0:np)
DOUBLE PRECISION PI
INTEGER NEXAMP
COMMON /NEXAMP /NEXAMP
PI = DACOS(-1.0D0)
C
C......................................................................
C
IF (NEXAMP .eq. 1) GOTO 1
IF (NEXAMP .eq. 2) GOTO 2
IF (NEXAMP .eq. 3) GOTO 3
IF (NEXAMP .eq. 4) GOTO 4
IF (NEXAMP .eq. 5) GOTO 5
IF (NEXAMP .eq. 6) GOTO 6
C
C......................................................................
C
C NON OSCILLATING SOLUTION
C
1 u = DSIN(PI*tx) * DSIN(PI*ty) * 1000.0D0
GOTO 10
C
C......................................................................
C
C HIGHLY OSCILLATING SOLUTION
C
2 u = DSIN(20.0d0*PI*tx) * DSIN(20.0d0*PI*ty)
GOTO 10
C
C......................................................................
C
C SOLUTION INCREASING SLOWLY IN TIME
C
3 u = Dexp(tx+ty) * 100000.0d0
GOTO 10
C
C......................................................................
C
C SOLUTION INCREASING RAPIDLY IN TIME
C
4 u = Dexp(3.0d0*(tx+ty))
GOTO 10
C
C......................................................................
C
C POLYNOMIAL IN X AND Y AS SOLUTION
C
5 u = tx**2 + ty**2
GOTO 10
C
C......................................................................
C
C ZERO SOLUTION
C
6 u = 0.0d0
GOTO 10
C
C......................................................................
C
10 CONTINUE
END
C
C
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C
C
subroutine F2(u,tx,ty,np)
C
C**********************************************************************
C
C INITIAL CONDITION DU/DT(X,Y,0) = F2(X,Y)
C
C**********************************************************************
C
integer np
DOUBLE PRECISION u(0:np,0:np), tx(0:np,0:np), ty(0:np,0:np)
DOUBLE PRECISION PI
INTEGER NEXAMP
COMMON /NEXAMP /NEXAMP
PI = DACOS(-1.0D0)
C
C......................................................................
C
IF (NEXAMP .eq. 1) GOTO 1
IF (NEXAMP .eq. 2) GOTO 2
IF (NEXAMP .eq. 3) GOTO 3
IF (NEXAMP .eq. 4) GOTO 4
IF (NEXAMP .eq. 5) GOTO 5
IF (NEXAMP .eq. 6) GOTO 6
C
C......................................................................
C
C NON OSCILLATING SOLUTION
C
1 u = 0.0d0
GOTO 10
C
C......................................................................
C
C HIGHLY OSCILLATING SOLUTION
C
2 u = 0.0d0
GOTO 10
C
C......................................................................
C
C SOLUTION INCREASING SLOWLY IN TIME
C
3 u = dsqrt(2.0d0) * dexp(tx+ty) * 100000.0d0
GOTO 10
C
C......................................................................
C
C SOLUTION INCREASING RAPIDLY IN TIME
C
4 u = 3.0d0 * dsqrt(2.0d0) * dexp(3.0d0*(tx+ty))
GOTO 10
C
C......................................................................
C
C POLYNOMIAL IN X AND Y AS SOLUTION
C
5 u = 0.0d0
GOTO 10
C
C......................................................................
C
C ZERO SOLUTION
C
6 u = 0.0d0
GOTO 10
C
C......................................................................
C
10 CONTINUE
END
C
C
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C
C
subroutine G(u,tx,ty,np,t)
C
C**********************************************************************
C
C BOUNDARY CONDITION U(X,Y,T) = G(X,Y,T)
C
C**********************************************************************
C
integer np
DOUBLE PRECISION u(0:np,0:np), tx(0:np,0:np), ty(0:np,0:np)
DOUBLE PRECISION t
DOUBLE PRECISION PI
INTEGER NEXAMP
COMMON /NEXAMP /NEXAMP
PI = DACOS(-1.0D0)
C
C......................................................................
C
IF (NEXAMP .eq. 1) GOTO 1
IF (NEXAMP .eq. 2) GOTO 2
IF (NEXAMP .eq. 3) GOTO 3
IF (NEXAMP .eq. 4) GOTO 4
IF (NEXAMP .eq. 5) GOTO 5
IF (NEXAMP .eq. 6) GOTO 6
C
C......................................................................
C
C NON OSCILLATING SOLUTION
C
1 u(0,0:np) = dsin(pi*tx(0,0:np)) * dsin(pi*ty(0,0:np))
+ * dcos(dsqrt(2.0d0)*pi*t) * 1000.0d0
u(np,0:np) = dsin(pi*tx(np,0:np)) * dsin(pi*ty(np,0:np))
+ * dcos(dsqrt(2.0d0)*pi*t) * 1000.0d0
u(1:np-1,0) = dsin(pi*tx(1:np-1,0)) * dsin(pi*ty(1:np-1,0))
+ * dcos(dsqrt(2.0d0)*pi*t) * 1000.0d0
u(1:np-1,np) = dsin(pi*tx(1:np-1,np)) * dsin(pi*ty(1:np-1,np))
+ * dcos(dsqrt(2.0d0)*pi*t) * 1000.0d0
GOTO 10
C
C......................................................................
C
C HIGHLY OSCILLATING SOLUTION
C
2 u(0,0:np) = dsin(20.0d0*pi*tx(0,0:np))
+ * dsin(20.0d0*pi*ty(0,0:np))
+ * dcos(20.0d0*dsqrt(2.0d0)*pi*t)
u(np,0:np) = dsin(20.0d0*pi*tx(np,0:np))
+ * dsin(20.0d0*pi*ty(np,0:np))
+ * dcos(20.0d0*dsqrt(2.0d0)*pi*t)
u(1:np-1,0) = dsin(20.0d0*pi*tx(1:np-1,0))
+ * dsin(20.0d0*pi*ty(1:np-1,0))
+ * dcos(20.0d0*dsqrt(2.0d0)*pi*t)
u(1:np-1,np) = dsin(20.0d0*pi*tx(1:np-1,np))
+ * dsin(20.0d0*pi*ty(1:np-1,np))
+ * dcos(20.0d0*dsqrt(2.0d0)*pi*t)
GOTO 10
C
C......................................................................
C
C SOLUTION INCREASING SLOWLY IN TIME
C
3 u(0,0:np) = dexp(tx(0,0:np)+ty(0,0:np))
+ * dexp(dsqrt(2.0d0)*t) * 100000.0d0
u(np,0:np) = dexp(tx(np,0:np)+ty(np,0:np))
+ * dexp(dsqrt(2.0d0)*t) * 100000.0d0
u(1:np-1,0) = dexp(tx(1:np-1,0)+ty(1:np-1,0))
+ * dexp(dsqrt(2.0d0)*t) * 100000.0d0
u(1:np-1,np) = dexp(tx(1:np-1,np)+ty(1:np-1,np))
+ * dexp(dsqrt(2.0d0)*t) * 100000.0d0
GOTO 10
C
C......................................................................
C
C SOLUTION INCREASING RAPIDLY IN TIME
C
4 u(0,0:np) = dexp(3.0d0* (tx(0,0:np)+ty(0,0:np)) )
+ * dexp(3.0d0*dsqrt(2.0d0)*t)
u(np,0:np) = dexp(3.0d0* (tx(np,0:np)+ty(np,0:np)) )
+ * dexp(3.0d0*dsqrt(2.0d0)*t)
u(1:np-1,0) = dexp(3.0d0* (tx(1:np-1,0)+ty(1:np-1,0)) )
+ * dexp(3.0d0*dsqrt(2.0d0)*t)
u(1:np-1,np) = dexp(3.0d0* (tx(1:np-1,np)+ty(1:np-1,np)) )
+ * dexp(3.0d0*dsqrt(2.0d0)*t)
GOTO 10
C
C......................................................................
C
C POLYNOMIAL IN X AND Y AS SOLUTION
C
5 u(0,0:np) = tx(0,0:np)**2 + ty(0,0:np)**2 + 2.0d0*t**2
u(np,0:np) = tx(np,0:np)**2 + ty(np,0:np)**2 + 2.0d0*t**2
u(1:np-1,0) = tx(1:np-1,0)**2 + ty(1:np-1,0)**2 + 2.0d0*t**2
u(1:np-1,np) = tx(1:np-1,np)**2 + ty(1:np-1,np)**2 + 2.0d0*t**2
GOTO 10
C
C......................................................................
C
C ZERO SOLUTION
C
6 u(0,0:np) = 0.0d0
u(np,0:np) = 0.0d0
u(1:np-1,0) = 0.0d0
u(1:np-1,np) = 0.0d0
C
C......................................................................
C
10 CONTINUE
END